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ABSTRACT 


We  report  progress  on  the  development  of  methods  in  a  number  of  specific  areas  of 
application  including  static,  non-cooperative  games  related  to  counter-  and  counter-counter- 
electromagnetic  interrogation  of  targets,  modeling  of  complex  viscoelastic  polymeric  mate¬ 
rials,  stochastic  and  deterministic  models  for  complex  networks  and  development  of  inverse 
problem  methodologies  (generalized  sensitivity  functions;  asymptotic  standard  errors)  for 
estimation  of  infinite  dimensional  functional  parameters  including  probability  measures  and 
temporal/spatial  dependent  functions  in  complex  nonlinear  dynamical  systems.  These  ef¬ 
forts  are  part  of  our  fundamental  research  in  a  modeling,  estimation  and  control  methodology 
(theoretical,  statistical  and  computational)  for  systems  in  the  presence  of  major  model  and 
observation  uncertainties. 
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1  Brief  Overview  Summary 

One  goal  of  a  principle  component  of  our  research  has  been  multiscale  issues  in  complex 
materials  as  a  means  to  understand  the  macroscopic  mechanical  properties  of  filled  rub¬ 
ber/elastomers  and  polymers  (including  biotissue)  and  their  connection  to  the  microscopic 
origin  of  stress  in  these  materials.  During  the  period  of  this  grant  we  have  made  significant 
progress  as  detailed  in  the  publications  [P13,P15].  In  these  efforts  we  incorporated  a  Rouse 
model  for  polymer  chains  into  the  linear  continuous  stick-slip  molecular-based  tube  repta- 
tion  ideas  of  Doi-Edwards  and  Johnson-Stacer.  This  treats  the  physically  constrained  (PC) 
molecular  stretches  as  internal  strain  variables  for  the  overall  PC/chemically  cross-linked 
(CC)  system.  It  yields  an  explicit  system  of  stress-strain  equations  for  the  system  permit¬ 
ting  simple  calculations  of  complex  stress-strain  relations.  The  model  that  as  developed  in 
[P13,P15]  treats  the  PC  molecule  as  entrapped  within  a  constraining  tube,  which  is  com¬ 
prised  of  both  CC  and  PC  molecules.  The  model  was  tested  and  validated  with  experimental 
data  sets  from  the  literature  including  data  from  polyisoprene,  and  polyisoprene  with  carbon 
black  reinforcement.  The  ideas  were  applied  in  [P 14]  to  model  shear  waves  in  biotissue  as 
part  development  of  inverse  problem  methodology  in  noninvasive  medical  diagnostics. 

We  continue  to  make  significant  progress  in  our  work  on  development  of  a  probability 
based  framework  (the  Prohorov  Metric  Framework  or  PMF)  to  treat  uncertainty  in  dy¬ 
namical  models.  Our  work  on  related  sensitivity  concepts  for  dynamical  systems  depending 
explicitly  on  the  probability  distributions  characterizing  modeling  heterogeneities  and  uncer¬ 
tainties  has  resulted  in  new  and  fundamental  theories  in  [P1,P7,P12].  Because  these  generally 
unknown  ’’parameters”  (the  probability  distributions  or  measures)  do  not  belong  to  a  Ba¬ 
nach  space,  fundamental  questions  regarding  sensitivity  derivatives  arise.  Our  progress  on 
theoretical  and  computational  fronts  by  was  achieved  by  treating  the  distributions  as  ele¬ 
ments  of  a  convex  subset  of  a  metric  space  with  the  Prohorov  metric;  sensitivities  are  defined 
in  terms  of  directional  derivatives.  We  demonstrated  the  power  of  these  ideas  on  both  linear 
and  nonlinear  systems  with  distribution  dependent  dynamics  as  outlined  in  [P7] .  Important 
applications  (in  modeling,  estimation  and  control  problems)  include  modeling  uncertainties 
in  complex  networks  such  as  food  production  [P16]  ,  communication,  logistical  supply  chains, 
etc.;  electromagnetic  systems  in  complex  dielectrics  [P2,P5,P9];  and  viscoelastic  materials 
[P13,P15]  as  well  as  models  for  production  of  biological  counter  measures  to  viral  and  toxic 
attacks  on  populations  [P7].  Theoretical  and  computational  efforts  based  on  [P4]  for  an 
asymptotic  distribution  theory  (for  large  sample  sizes)  are  currently  being  pursued  with 
slow  but  steady  progress.  Some  of  the  theoretical  efforts  are  motivated  by  our  continuing 
efforts  to  understand  dispersal  of  pulsed  electromagnetic  energy  in  interrogation  problems 
[P3,P10,P11,P18,P19]  for  complex  materials  and  targets.  Our  efforts  on  understanding  the 
use  of  standard  statistical  methodologies  in  finite  dimensional  inverse  problems  [P6,P8,P  1 7] 
is  of  significant  importance  to  this  pursuit  of  new  uncertainty  methodologies. 

Finally,  we  are  able  report  major  progress  in  our  long  term  efforts  on  anti-interrogation 
and  anti-anti-interrogation  methodologies  using  controllable  ferroelectric  and  ferromagnetic 
layers  coatings  on  conducting  objects  such  as  airfoils  and  missiles  to  provide  an  attenuation 
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capability  against  electromagnetic  interrogation.  Our  early  computational  findings  [P3]  led 
to  investigation  of  a  new  class  of  problems  wherein  problems  are  formulated  as  a  zero-sum  two 
player  minmax  game.  Here  each  player  (evader  and  interrogator)  must  choose  a  best  strategy 
in  the  presence  of  uncertainty  as  to  what  the  other  player  will  do.  In  [P5]  we  successfully 
formulated  static  minmax  games  over  spaces  of  probability  measures  where  the  full  power 
of  the  Prohorov  Metric  Framework  and  our  associated  theory,  approximation,  sensitivity 
and  computational  methodology  plays  an  essential  enabling  role.  We  have  now  turned  to 
dynamic  games  where  both  the  evader  and  interrogator  has  some  real  time  dynamic  control 
over  his  own  uncertain  strategies  but  has  only  incomplete  information  about  his  adversary’s 
strategies;  these  are  modeled  as  Markov  diffusion  processes  at  the  moment.  Our  team  is 
making  great  progress  in  this  direction. 

1.1  Publications  supported  in  part  under  this  grant 

[PI]  H.T.  Banks  and  H.K.  Nguyen,  Sensitivity  of  dynamical  systems  to  Banach  space  pa¬ 
rameters,  J.  Math.  Analysis  and  Applications,  323  (2006),  141  161. 

[P2]  H.T.  Banks  and  N.L.  Gibson,  Electromagnetic  inverse  problems  involving  distribu¬ 
tions  of  dielectric  mechanisms  and  parameters,  Quarterly  Appl.  Math.,  64  (2006),  749  795. 

[P3]  H.T.  Banks,  K.  Ito  and  J.A.  Toivanen,  Determination  of  interrogating  frequencies  to 
maximize  electromagnetic  backscatter  from  objects  with  material  coatings,  Communications 
in  Comp.  Physics,  1  (2006),  357-377. 

[P4]  H.T.  Banks  and  J.L.  Davis,  A  comparison  of  approximation  methods  for  the  esti¬ 
mation  of  probability  distributions  on  parameters,  CRSC-TR05-38,  October,  2005;  Applied 
Numerical  Mathematics,  57  (2007),  753  777. 

[P5]  H.T.  Banks,  S.L.  Grove,  K.  Ito  and  J.A.  Toivanen,  Static  two-player  evasion- 
interrogation  games  with  uncertainty,  CRSC-TR06-16,  June,  2006;  Comp,  and  Applied 
Math,  to  appear. 

[P6  ]  H.T.  Banks,  Stacey  L.  Ernstberger  and  Sarah  L.  Grove,  Standard  errors  and  con¬ 
fidence  intervals  in  inverse  problems:  Sensitivity  and  associated  pitfalls,  J.  Inverse  and  Ill- 
posed  Problems,  15  (2007),  1-18. 

[P7]  H.T.  Banks,  S.  Dediu  and  H.K. Nguyen,  Time  delay  systems  with  distribution  de¬ 
pendent  dynamics,  IFAC  Annual  Reviews  in  Control,  31  (2007),  17  26. 

[P8]  H.T.  Banks,  G.M.  Kepler,  H.K.  Nguyen  and  J.  Webster-Cyriaque),  Inverse  problems 
and  model  validation:  An  example  from  latent  virus  reactivation  CRSC-TR06-18,  August, 
2006;  J.  Inverse  and  Ill-posed  Problems,  15  (2007),  19-35. 

[P9  ]  H.T.  Banks,  V.A.  Bokil  and  N.L.  Gibson,  Parameter  estimation  versus  homogeniza¬ 
tion  techniques  in  time-domain  characterization  of  composite  dielectrics,  CRSC-TR06-20, 
August,  2006;  J.  Inverse  and  Ill-posed  Problems,  15  (2007),  117  135. 

[P10  ]  H.T.  Banks,  V.A.  Bokil  and  N.L.  Gibson,  Analysis  of  stability  and  dispersion  in 
a  finite  element  method  for  Debye  and  Lorentz  dispersive  media,  CRSC-TR06-21,  August, 
2006;  Numerical  Methods  for  Partial  Differential  Equations,  submitted. 
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TR07-08,  March,  2007;  IEEE  Transactions  on  Biomedical  Engineering,  submitted. 

[P18]  K.  Ito,  Q.  Huynh  and  J.  Toivanen,  A  fast  Helmholtz  solver  for  scattering  by  a 
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[P19]  K.Ito  and  J.  Toivanen,  A  fast  iterative  solver  for  scattering  by  elastic  objects  in 
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•  V.A.  Bokil,  Post  doc,  North  Carolina  State  University  (now  at  Oregon  State  U.) 

•  S.  Dediu,  Post  doc,  North  Carolina  State  University 

•  S.  Hu,  Post  doc,  North  Carolina  State  University  N.  Luke,  Post  doc,  North  Carolina 
State  University  (now  at  EPA) 

•  H.  Nguyen,  Post  doc,  North  Carolina  State  University  (now  at  Ctr  for  Naval  Analysis) 

•  S.N.  Wynne,  Post  doc,  North  Carolina  State  University 

•  J.  Ernstberger,  Graduate  Student  at  North  Carolina  State  University 

•  S.  Grove,  Graduate  Student  at  North  Carolina  State  University 

•  J.  Hood,  Graduate  Student  at  North  Carolina  State  University 

•  S.  Joyner,  Graduate  Student  at  North  Carolina  State  University 

•  N.  Luke,  Graduate  Student  at  North  Carolina  State  University 

•  J.  Samuels,  Graduate  Student  at  North  Carolina  State  University 

•  L.  Stroud/Dick,  Graduate  Student  at  North  Carolina  State  University 
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•  K.  Sutton,  Graduate  Student  at  N  C  State  University  and  Arizona  State  U. 

•  W.  Wan,  Graduate  Student  at  North  Carolina  State  University 

•  C.  Hong,  Graduate  Student  at  North  Carolina  State  University 

3  Relevance  of  Research  to  DOD/AF  Missions 

The  development  of  dynamic  electromagnetic  pursuit/evasion  games  with  uncertainty  and/or 
stochasticity  is  fundamental  to  DOD/AF  development  of  an  electromagnetic/mechanical 
active  stealth  technology  for  aerospace,  land  and  marine  vehicles  and  other  potential  targets 
of  electromagnetic  imaging.  The  origins  of  this  proposed  effort  lies  in  long  term  collaborative 
efforts  with  Dr.  R.  Albanese  and  his  group  at  AFRL,  Brooks  City  Base  in  San  Antonio. 
Indeed  many  of  the  issues  we  address  are  in  direct  response  to  questions  raised  by  Dr. 
Albanese  and  his  team.  We  have  enjoyed  a  continued  close  collaborative  relationship  with 
the  group  during  the  efforts. 

The  development  of  inverse  problem  methodology  for  uncertainty  quantification  and  the 
development  of  general  classes  of  nonlinear  complex  nodal  network  models  with  inherent  un¬ 
certainties  is  directly  relevant  to  several  recently  announced  “Discovery  Challenge  Thrusts” 
in  AFOSR-BAA-2007-08.  In  particular,  our  efforts  are  related  to  the  programs  announced 
on  “  Robust  Decision  Making”  and  “Complex  Networked  Systems”  which  are  concerned  with 
operational  environments  where  “unexpected  events,  attacks  on  and  degradations  of  the  net¬ 
work  centric  system  and  sudden  changes  in  goals  and  objectives  ”  must  be  accommodated. 
The  need  for  new  modeling  paradigms  which  include  robustness  in  estimation  and  control 
in  models  with  uncertainty  is  well  recognized. 

We  expect  our  previous  efforts  to  be  of  direct  interest  to  scientists  and  engineers  in  the 
Information  Directorate  at  AFRL,  Rome,  NY.  Indeed,  this  group  has  recently  hired  one 
of  Banks’  Ph.  D.  students  and  co-authors,  Sarah  Grove,  whose  thesis  was  completed  in 
September,  2007  and  who  worked  on  the  electromagnetic  pursuit/evasion  games.  We  hope 
to  initiate  collaboration  with  their  investigators,  especially  in  the  nodal  network  modeling 
aspects  of  our  efforts. 


4  Detailed  Research  Summaries 

4.1  Interrogation/ Counter-Interrogation/ 

Counter-Counter-Interrogation  (I/CI/CCI) 

In  continuing  efforts  [17,  21,  22]  in  collaborative  efforts  with  Dr.  R.A.  Albanese  (AFRL, 
Brooks  AFB)  and  his  associates,  we  have  made  significant  progress  on  counter-interrogation 


7 


methodologies.  We  have  assumed  controllable  ferroelectric  and  ferromagnetic  layers  coating 
a  conducting  object  to  provide  an  attenuation  capability  against  electromagnetic  interroga¬ 
tion.  The  problems  can  be  formulated  as  robust  optimization  and/or  two-player  games.  We 
have  shown  that  the  scattered  field  due  to  interrogation  can  be  attenuated  even  with  the 
assumption  of  uncertainty  in  the  interrogation  wave  numbers.  The  controllable  layer  com¬ 
posed  of  ferromagnetic  and  ferroelectric  materials  can  be  incorporated  into  a  mathematical 
formulation  based  on  the  time-harmonic  Maxwell  equation.  Fresnel’s  law  for  the  reflectance 
index  has  been  extended  to  the  electromagnetic  propagation  in  anisotropic  composite  lay¬ 
ers  of  ferromagnetic  and  electronic  devices  and  used  to  demonstrate  feasibility  of  control 
of  reflections.  Our  methodology  has  also  been  tested  for  a  non-planar  geometries  of  the 
conducting  object  (an  NACA  airfoil  and  a  missile)  in  which  we  reported  our  findings  in  the 
form  of  reduced  radar  cross  sections  (RCS). 

We  subsequently  [22]  have  shown  that  even  rather  simple  counter-counter-interrogation 
measures  by  the  interrogator  will  defeat  straight  forward  evasion  design.  Thus  we  were  led 
to  investigations  of  how  each  of  the  evader  and  the  interrogator  might  try  to  confuse  the 
other  by  adding  uncertainty  in  their  design  strategies.  Our  early  computational  findings  led 
to  investigation  of  a  new  class  of  problems  wherein  problems  are  formulated  as  a  zero-sum 
two  player  minmax  game.  Here  each  player  must  choose  a  best  strategy  in  the  presence  of 
uncertainty  as  to  what  the  other  player  will  do.  This  leads  in  [17]  to  minmax  games  over 
spaces  of  probability  measures  where  the  full  power  of  the  Prohorov  Metric  Framework  and 
our  associated  theory,  approximation,  sensitivity  and  computational  methodology  plays  an 
essential  enabling  role. 

Specifically,  in  [21]  we  demonstrated  that  it  is  possible  to  design  ferroelectric  materials 
with  appropriate  dielectric  permittivity  and  magnetic  permeability  to  significantly  attenuate 
reflections  of  electromagnetic  interrogation  signals  from  highly  conductive  targets  such  as 
airfoils  and  missiles.  This  was  done  under  assumptions  that  the  interrogating  input  signal 
is  uniformly  likely  to  come  from  a  sector  of  interrogating  angles  a  €  [ao,  au]  (a  =  |  —  4> 
where  4>  is  the  angle  of  incidence  of  the  input  signal)  but  that  the  evader  has  knowledge 
of  the  interrogator’s  input  frequency  or  frequencies  (denoted  here  as  the  interrogator  design 
frequencies  Id)-  Example  radar  cross  sections  defined  by 

RCS  (a)  =  101og10  ^|FQ(a  +  7r)|2^ 

are  shown  for  two  optimized  coating  materials  vs.  that  for  no  coating  in  Figure  1.  (Here 
Fa(a+ 7r)  is  the  reflected  far  field  for  incident  angle  a  based  on  inverse  scattering  for  Maxwell’s 
equations  as  explained  in  detail  in  [21].)  Corresponding  reflected  field  intensities  for  compar¬ 
ison  between  the  no  coating  layer  case  and  the  optimized  complex  parameter  case  of  Figure 
1  are  depicted  in  Figure  2  for  an  angle  of  interrogation  a  =  n/4. 
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270 


Figure  1:  The  RCS  for  optimized  complex-valued  material  parameters  (e*,/z*)  (blue  line); 
for  optimized  dielectric  only  (e*,/xr  =  1)  (green  line);  and  for  no  coating  (er  =  1;  [iT  =  1) 
(red  line)  for  wavelength  A  =  1/4. 


Figure  2:  The  reflected  field  E ^  for  no  coating  [tr  =  l,/xr  =  1)  (left)  and  for  the  case 
of  optimized  complex-valued  material  parameters  (e*,/u*)  (right)  for  A  =  1/4  and  angle  of 
interrogation  a  =  7r/4. 
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These  results  were  further  sharpened  and  illustrated  in  [22]  where  a  series  of  different 
material  designs  were  considered  to  minimize  over  a  given  set  of  input  design  frequencies  Ip 
the  maximum  reflected  field  from  input  signals.  In  addition,  a  second  critical  finding  was 
obtained  in  that  it  was  shown  that  if  the  evader  employed  a  simple  counter  interrogation 
design  based  on  a  fixed  set  (assumed  known)  of  interrogating  frequencies  Ip,  then  by  a 
rather  simple  counter-counter  interrogation  strategy  (use  an  interrogating  frequency  little 
more  than  10%  different  from  the  assumed  design  frequencies),  the  interrogator  can  easily 
defeat  the  evader’s  material  coatings  counter  interrogation  strategy  to  obtain  strong  reflected 
signals. 

From  the  combined  results  of  [21,  22]  it  is  thus  rather  easily  concluded  that  the  evader 
and  the  interrogator  must  each  try  to  confuse  the  other  by  introducing  significant  uncertainty 
in  their  design  and  interrogating  strategies,  respectively.  This  concept,  which  we  refer  to 
as  mixed  strategies  in  recognition  of  previous  contributions  to  the  literature  on  games  (von 
Neumann’s  finite  mixed  strategies- see[17]),  leads  to  two  player  non-cooperative  games  with 
probabilistic  strategy  formulations.  These  can  be  mathematically  formulated  as  two  sided 
optimization  problems  over  spaces  of  probability  measures,  i.e.,  minmax  games  over  sets  of 
probability  measures.  In  [17]  we  provided  a  mathematically  precise  formulation  of  such  a 
class  of  two  sided  optimization  problems  and  discussed  our  initial  computational  efforts  on 
such  problems.  Approximation  methods  were  introduced  and  their  computational  efficacy 
were  demonstrated  with  several  simple  examples. 

More  precisely,  we  considered  electromagnetic  interrogation  of  objects  in  the  context  of 
minmax  evader-interrogator  games  where  each  player  has  uncertain  information  about  the 
adversary’s  capabilities.  The  minmax  cost  functional  is  based  on  reflected  fields  from  an 
object  such  as  an  airfoil  or  missile.  We  used  the  simplest  reflection  coefficient  based  on  a 
simple  planar  geometry  (e.g.,  see  Figure  3)  using  Fresnel’s  formula  for  a  perfectly  conducting 
half  plane  which  has  a  coating  layer  of  thickness  d  with  dielectric  permittivity  e  and  magnetic 
permeability  p.  A  normally  incident  (<p  =  0)  electromagnetic  wave  with  the  frequency  /  is 
assumed  to  impinge  the  half  plane.  Then  the  corresponding  wavelength  A  in  air  is  A  =  c/ /, 
where  the  speed  of  light  is  c  =  0.3  x  109. 

The  reflection  coefficient  R  for  the  wave  is  given  by 


R  = 


a  +  b 
1  +  ab' 


(4.1) 


where 

a=i-v^  ^  b  =  e4in^fifd/c  (4.2) 

e  +  y/ep 

This  expression  can  be  derived  directly  from  Maxwell’s  equation  by  considering  the  ratio  of 
reflected  to  incident  wave  for  example  in  the  case  of  parallel  polarized  ( TEX )  incident  wave 
(see  [21,  54]).  An  alternative  and  much  more  computationally  intensive  approach  (which 
may  be  necessitated  by  some  target  geometries)  employs  the  far  field  pattern  (see  [38])  for 
intensity  BR(e,  p,  a,  f)  of  reflected  waves  computed  directly  using  Maxwell’s  equations  (see 
[21,  22])  and  allowing  a  non-normal  angle  of  incidence  a  =  |  —  0  ^ 
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Ambient 


Figure  3:  Interrogating  high  frequency  wave  impinging  (angle  of  incidence  0)  on  coated 
(thickness  d)  perfectly  conducting,  surface 

The  evader  and  the  interrogator  are  each  assumed  subject  to  uncertainties  as  to  the 
actions  of  the  other.  The  evader  wants  to  choose  a  best  coating  design  (i.e.,  best  e’s  and 
p's  )  while  the  interrogator  wants  to  choose  best  angles  of  interrogation  a  and  interrogating 
frequencies  /.  Each  player  must  act  in  the  presence  of  incomplete  information  about  the 
other’s  action.  Partial  information  regarding  capabilities  and  tendencies  of  the  adversary 
can  be  embodied  in  probability  distributions  for  the  choices  to  be  made.  That  is,  we  may 
formalize  this  by  assuming  the  evader  may  choose  (with  an  as  yet  to  be  determined  set 
of  probabilities)  dielectric  permittivity  and  magnetic  permeability  parameters  (e,p)  from 
admissible  sets  £xA4  while  the  interrogator  chooses  angles  of  interrogation  and  interrogating 
frequencies  (a,  /)  from  sets  A  x  T .  The  formulation  here  is  based  on  the  mixed  strategies 
proposals  of  von  Neumann  [4,  60,  61]  and  the  ideas  can  be  summarized  as  follows.  The 
evader  does  not  choose  a  single  coating,  but  rather  has  a  set  of  possibilities  available  for 
choice.  He  only  chooses  the  probabilities  with  which  he  will  employ  the  materials  on  a 
target.  This,  in  effect,  disguises  his  intentions  from  his  adversary.  By  choosing  his  coatings 
randomly  (according  to  a  best  strategy  to  be  determined  in,  for  example,  a  minmax  game), 
he  prevents  adversaries  from  discovering  which  coating  he  will  use  indeed,  even  he  does  not 
know  which  coating  will  be  chosen  for  a  given  target.  The  interrogator,  in  a  similar  approach, 
determines  best  probabilities  for  choices  of  frequency  and  angle  in  the  interrogating  signals. 
Note  that  such  a  formulation  tacitly  assumes  that  the  adversarial  relationship  persists  with 
multiple  attempts  at  evasion  and  detection. 

The  associated  minmax  problem  consists  of  the  evader  choosing  distributions  Pe(e,  p)  over 
£xMto  minimize  the  reflected  field  while  the  interrogator  chooses  distributions  P,(a,  /)  over 
Ax  V  to  maximize  the  reflected  field.  If  BR(e,  p,  a,  f)  is  the  chosen  measure  of  reflected 
field  and  Ve  =  V(£  x  A4)  and  Vi  =  V(A  x  V)  are  the  corresponding  sets  of  probability 
distributions  or  measures  over  £  x  Ad  and  A  x  V,  respectively,  then  the  cost  functional  for 
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the  minmax  problem  can  be  defined  by 

J(Pe,  Pi)  =  f  [  |  BR(e,  ijl,  a,  f)\2d,Pe(e,  p)dP(a,  /).  (4.3) 

JSxM  JAxF 

The  problems  thus  formulated  are  special  cases  of  classical  static  zero  sum  two  player 
non-cooperative  games  [4,  33]  where  the  evader  minimizes  over  Pe  G  Ve  and  the  interrogator 
maximizes  over  Pi  €  Vi.  In  such  games  one  defines  upper  and  lower  values  for  the  game  by 

J  =  inf  sup  J(Pe,  Pi)  (4.4) 

P'ZPe  p.sp, 


and 


J  —  sup  inf  J(Pe,Pi). 

~  Pl€V, 


(4.5) 


The  first  represents  a  security  level  (worst  case  scenario)  for  the  evader  while  the  latter  is 
a  security  level  for  the  interrogator.  It  is  readily  argued  that  J  <  J  and  if  the  equality 
J*  =  J  =  J  holds,  then  J*  is  called  the  value  of  the  game.  Moreover,  if  there  exist  P*  €  Ve 
and  P*  €  Vi  such  that 


r  =  J(p;,  P*)  =  minP'€Ve  J(Pe ,  P*)  =  max  J(Pe*,  P4), 


then  (P*,  P*)  is  a  saddle  point  solution  or  non-cooperative  equilibrium  of  the  game. 

To  investigate  theoretical,  computational  and  approximation  issues  for  these  problems,  it 
is  necessary  to  put  a  topology  on  the  space  of  probability  measures.  Given  a  set  7,  a  natural 
choice  for  V(I)  is  the  Prohorov  metric  topology  [35,  51,  57]  as  used  for  minimization/inverse 
problems  in  [7,  9,  16,  30].  Prohorov  metric  (p)  convergence  is  weak*  convergence  in  V  when 
V  is  considered  as  a  subset  of  the  topological  dual  C*B{I)  of  the  space  Cb(I)  of  bounded 
continuous  functions  on  7.  More  precisely,  convergence  p(Pt,  P)  — ♦  0  in  the  Prohorov  metric 
is  equivalent  to 

J /(^)dPfe(O  >  JmdPiO  for  all  bounded,  uniformly  continuous  /  :  7  — >  R1. 

It  is  known  [7,  9]  that  if  7  is  a  complete  metric  space,  then  V  taken  with  the  Prohorov 
metric  is  a  complete  metric  space.  Moreover,  if  7  is  compact,  then  so  is  V.  Using  these 
properties  and  arguments  similar  to  those  in  [7,  15,  16,  24,  30],  one  can  develop  well-posedness 
and  approximation  results  for  the  minmax  problems  defined  above.  Efficient  computational 
methods  that  correspond  to  von  Neumann’s  finite  mixed  strategies  [4]  can  readily  be  presented 
in  this  context.  These  can  be  based  on  several  approximation  theories  that  have  been 
recently  developed  and  used.  The  first,  developed  in  [7]  and  based  on  Dirac  delta  measures, 
is  essentially  a  zero-order  spline  or  Walsh  element  scheme,  while  the  second,  developed  and 
used  in  [30],  is  a  first-order  or  piecewise  linear  spline  or  finite  element  scheme.  These  are 
discussed  in  more  detail  in  [17]. 
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To  establish  existence  of  a  saddle  point  solution  for  the  evasion-interrogation  problems 
formulated  above,  one  can  employ  a  fundamental  result  of  von  Neumann  [60,  61]  as  stated 
by  Aubin  ([4]-see  p.  126).  As  detailed  in  [17],  one  can  argue  a  desired  well-posedness  results. 

Suppose  £,  Ai,  4>,  T  are  compact  (a  reasonable  assumption  for  the  applications  of  interest 
to  us  here)  and  define  the  spaces  Xo  =  V(£  x  Ad),  Vo  =  V($  x  T)  taken  with  the  Prohorov 
metric.  Then  X0,  Y0  are  compact,  convex  subsets  of  X  =  C*B{£  x  M.)  and  Y  =  C*B($  x  T), 
respectively.  Moreover,  there  exists  (P*,  P*)  E  V{£  x  M)  x  V($  x  T)  such  that 


j(p;,pn 


min  max  J(Pe,Pi ) 

V(£xM)V(' PxX) 


max  min 
V(<t>xF)V(£xM) 


J (Pe,  Pi)- 


In  light  of  equations  (4.4)  and  (4.5),  this  important  results  allows  each  player  to  inves¬ 
tigate  worst  case  scenarios  (assuming  he  optimizes  his  own  situation)  even  in  presence  of 
uncertainty  about  his  adversary’s  action. 

In  [17]  we  also  developed  an  approximation  framework  and  a  corresponding  convergence 
theory  for  classes  of  min-max  problems  based  on  approximations  of  the  measures  using 
the  delta  measures  of  [7].  We  considered  a  number  of  different  material  designs  from  [22] 
that  were  developed  earlier  in  investigations  of  counter  interrogation  and  counter-counter 
interrogation  feasibility.  We  used  a  number  of  these  materials  to  define  families  of  admissible 
evader  sets  £  x  A4  in  our  initial  exploratory  calculations  in  demonstrating  ease  of  calculations 
for  problems  as  outlined  here.  Details  are  given  in  [17]. 


4.2  Inverse  Problems  for  Systems  with  Functional  Parameters 

In  [8]  we  considered  a  class  of  probability  measure  dependent  dynamical  systems  such  as 
those  that  arise  in  the  study  of  multiscale  phenomena  in  diverse  fields  such  as  immunological 
population  dynamics,  viscoelasticity  of  polymers  and  rubber,  and  polarization  in  dielectric 
materials.  We  established  conditions  for  existence  and  uniqueness  of  the  forward  problem 
and  well-posedness  (including  method  stability  under  numerical  approximations)  for  the 
inverse  problem  of  estimating  the  probability  measures.  Moreover,  we  presented  a  general 
theoretical  framework  including  implementable  approximation  ideas  for  inverse  problems 
involving  measure  dependent  dynamical  systems.  This  framework  dealt  with  systems  in 
which  the  expected  state  dynamics  are  given  by  an  abstract  differential  equation 

x(t)  =  F{t,x{t),P) , 

where  the  right  side  is  dependent  on  the  probability  measure  P.  In  general,  the  system  can 
also  represent  a  delay  or  partial  differential  equation. 

The  motivation  for  our  efforts  include  several  important  examples  of  measure  dependent 
dynamics  that  have  arisen  in  recent  applications  that  are  of  direct  interest  to  AFOSR,  AFRL, 
and  other  DOD  agencies.  Realistic  models  in  viscoelasticity  of  polymers  and  rubber  as  well  as 
polarization  in  dielectric  materials  share  certain  features  with  our  formulation.  Specifically, 
recent  studies  [28,  29,  18]  of  molecular-based  stick-slip  reptation  models  for  heterogeneous 
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(4.6) 


viscoelastic  polymer  chain  materials  involve  measure  dependent  systems  of  the  form 


d2u(t,x)  d 


where  u  is  the  tensile  displacement,  <j  is  the  measure  dependent  stress 
cr(t,x,P)  =  ge(e(t),e(t))  +  v  J  e1(t,x;T)dP(r) , 


(4.7) 


with  strain  e  =  and  “internal”  strain  t\  defined  by 

-sr( t ,  x;  t )  +  -ei (t,  x ;  r)  =  e(t,  x)/i(c(t,  x)) .  (4.8) 

ar  r 

Here  the  measure  P  accounts  for  the  distribution  of  relaxation  properties  in  the  heteroge¬ 
nous  long  chain  polymer  molecules.  In  [25,  26]  it  is  shown  that  similar  models  (with  different 
nonlinearities  in  (4.7),  (4.8))  are  important  if  one  replaces  Fung  kernels  [43]  with  equations 
for  distributed  molecular  mechanisms  in  describing  shear  response  in  biotissue.  In  another 
important  application  [15,  16],  the  systems  are  the  usual  Maxwell’s  equations  for  the  elec¬ 
tromagnetic  fields  E  and  H  in  a  heterogeneous  dielectric  and  are  given  by 

V  x  E  =  -&■  V  ■  D  =  0 
VxH  =  ^  +  J  V  ■  //  =  0 
D  —  erE  P , 


with  the  exception  that  now  the  probability  measure  ( P )  dependent  macroscopic  polarization 
V  is  given  by 

P{t)  =  [  \pi(t]Ti)  +P2(t]T2)]dP(T1,T2). 

In  this  case,  the  microscopic  orientational  (Debye)  polarization  p\(t\T\)  is  defined  by 

Pi  +  — Pi  =  ZE , 

Tl 

and  the  microscopic  electronic  (Lorentz)  polarization  P2(^;t2)  is  defined  by 


mp 2  +  cp2  +  kp2  =  iE , 

with  t2  =  ^r.  In  the  equations  for  both  p\  and  p2,  the  parameters  T\  and  r2  represent 
relaxation  parameters  which  may  vary  over  the  admissible  set  T\  x  T2  according  to  some 
unknown  but  sought  after  probability  measure  P  =  P{t\  ,  t2). 

In  each  of  these  examples,  one  seeks  to  characterize  the  material  behavior  to  perturba¬ 
tions  by  finding  a  measure  P*  that  provides  the  “best”  mathematical  system  response  when 
compared  to  observations  of  the  physical  system.  Thus,  inverse  problems  involving  complex 
nonlinear  systems  with  function  space  (probability  measures,  densities,  as  well  as  space  and 
time  dependent  coefficients)  “parameters”  to  be  estimated  are  important  in  many  modern 
applications. 

As  part  of  our  ongoing  study  of  estimation  of  function  space  parameters  in  complex  non¬ 
linear  systems,  we  have  also  developed  a  sensitivity  theory  [13,  14,  27]  for  systems  depending 
explicitly  on  functions  and  probability  measures. 
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4.3  Traditional  Sensitivity 

We  consider  first  a  parameter  estimation  problem  for  the  general  nonlinear  dynamical  system 


x(t)  =  g(t,x(t),6) 
x(t0)  =  x0, 


with  discrete  time  observations 


Vj  =  f(tj, 6)  +  0  =  Cx(tj,9)  +  ej,  j  =  1, . . .  ,n  (4.10) 

where  RN ,  /,  ej  €  RM  and  6  G  Rp.  The  matrix  C  is  an  M  x  N  matrix  which  gives  the 
observation  data  in  terms  of  the  components  of  the  state  variable  x. 

Traditional  sensitivity  functions  (TSF)  are  classical  sensitivity  functions  used  in  mathe¬ 
matical  modeling  to  investigate  how  the  output  of  a  model  changes  when  the  parameters  and 
the  initial  conditions  vary.  In  order  to  quantify  the  variation  in  the  state  variable  x(t)  with 
respect  to  changes  in  the  parameter  6  and  the  initial  condition  xo,  we  are  naturally  led  to 
consider  the  (traditional)  sensitivity  functions  (TSF) 

k=l,...,p,  (4.11) 


1  =  1,. ..,N,  (4.12) 

where  x0/  is  the  Z-th  component  of  the  initial  condition  x0 .  When  the  function  g  is  sufficiently 
regular,  the  solution  x  is  differentiable  with  respect  to  9k  and  x0/,  and  therefore  our  sensitivity 
functions  sgk  and  rXQl  are  well  defined. 

Often  in  practice  the  model  under  investigation  is  simple  enough  to  allow  one  to  compute 
the  sensitivity  functions  (4.11)  and  (4.12)  directly.  However,  for  most  more  complex  models, 
it  is  often  preferable  (and  by  now  standard)  to  consider  these  sensitivity  functions  as  solutions 
of  their  own  dynamical  systems. 

From  the  sensitivity  analysis  theory  for  dynamical  systems  (numerous  references  are  given 
in  [13,  14,  27]),  we  have  that  s  =  (s^17 . . .  ,sep )  is  an  iV  x  p  vector  function  that  satisfies  the 
ODE  system 


and 


.  .  OT 

Sek^  =  Wk^' 


s(f)  =  gx(t,x,6)s(t)  +  gg{t,x,6),  (4.13) 

s(to)  =  o  jVxp- 

Similarly,  the  sensitivity  functions  with  respect  to  the  components  of  the  initial  condition  xq 
define  an  N  x  N  vector  function  r  =  (rI01, . . .  ,rl0N),  which  satisfies 

r(t)  =  gx(t,x,0)r(t),  (4.14) 

r(t0)  =  In*n- 
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Here  we  denote  by  gx  =  dg/dx  and  by  g$  =  dg/d6 ,  the  derivatives  of  g  with  respect  to  x 
and  6,  respectively. 

The  equations  (4.13)  and  (4.14)  are  used  in  conjunction  with  equation  (4.9)  to  numerically 
compute  the  sensitivities  s  and  r  for  general  cases  when  the  function  g  is  complicated  enough 
to  not  allow  a  closed  form  solution  by  direct  integration.  As  we  shall  explain  below,  these 
sensitivity  systems  are  important  components  of  a  standard  asymptotic  theory  (see  Chapter 
12  of  [58])  for  standard  errors  and  confidence  intervals. 

To  deal  with  sensitivity  for  complex  systems  with  general  function  space  parameters,  in 
[27]  we  considered  a  special  case  of  an  abstract  nonlinear  ODE  where  the  state  space  is  a 
general  Banach  space  X  and  the  parameter  space  M.c  is  a  convex  subset  of  a  Banach  space 

M.  such  as,  for  example,  M.  =  L2.  Specifically,  if  one  considers  a  general  nonlinear  ODE  of 

the  form 

*(t)  =  g(t,x(t),p),  (4.15) 

where  g  :  M+  x  X  x  Ad  — >  X  and  X  and  Ad  are  complex  Banach  spaces,  one  can  under 
reasonable  assumptions  establish  both  well-posedness  and  a  general  sensitivity  theory.  Let 
B(X ,  Y )  denote  the  space  of  bounded  linear  operators  from  X  onto  Y  and  suppose  the 
function  g(t,x,p)  of  (4.15)  has  continuous  Frechet  derivatives  gx(t,x,p)  with  respect  to  x 
and  gn(t,x,n)  with  respect  to  p.  Then  one  can  show  that  the  Frechet  derivative  y(t)  = 
T^x(t,to,xo,  g.)  exists  with  y(t)  in  B(Xi,X)  satisfying  the  equation 

y(t)  =  gx(t,x(t,t0,x0,p),p)y(t) +  g^(t,x(t,t0,x0,p),p), 
y(to)  =  0, 


for  t  >  to. 

With  the  parameter  space  of  probability  density  functions  Adc,  which  is  a  convex  subset 
of  a  Banach  space  Ad ,  this  sensitivity  theory  can  also  be  applied  using  directional  derivatives 
instead  of  the  Frechet  derivative.  However,  it  is  shown  in  [27]  that  the  directional  derivative 
of  a  continuous  function  g  is  the  Frechet  derivative  on  Ad  restricted  to  q  —  p  where  p,  q  €  Adc- 

In  order  to  accommodate  more  general  problems  with  Dirac  (discrete)  measures  or  mea¬ 
sures  with  a  continuous  component  and  a  saltus  component,  the  theoretical  results  above 
were  extended  to  a  general  convex  metric  space.  This  is  necessary  because  the  parameter 
space  associated  with  the  Prohorov  metric  is  no  longer  a  Banach  space  but  is  only  a  special 
case  of  a  convex  metric  space  (Adi,  Although  the  Prohorov  metric  is  not  conceptually 

easy  to  use,  it  generates  an  equivalent  topology  to  the  weak  L2  topology  (e.g.,  see  [30])  which 
is,  of  course,  the  same  as  the  weak*  topology  in  the  L2  case.  Therefore,  the  Prohorov  metric 
may  be  applied  in  numerical  approximation  in  distribution  dependent  problems  taking  ad¬ 
vantage  of  its  relation  in  convergence  to  the  weak  L 2  convergence.  To  treat  such  problems, 
in  [13]  we  extended  the  theoretical  results  mentioned  above  to  the  case  where  the  parameter 
space  is  a  convex  metric  space.  Let  (Adi,d>fi)  denote  a  general  convex  metric  space  with 
distance  dMl  and  X  denote  a  general  complex  Banach  space.  Again  we  consider  a  general 
nonlinear  abstract  ODE 


x(t)=  g(t,x{t),n  i),  x(0)  =  0, 


(4.16) 
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where  j:I+xXxMi  — >  X  is  continuous  in  all  three  variables  and  Frechet  differentiable  in 
x.  Here  the  solution  x  G  X  and  the  parameter  G  M.\.  The  conditions  for  and  statement  of 
existence  and  uniqueness  of  solutions  of  equation  (4.16)  along  with  continuous  dependence  of 
solutions  for  the  general  convex  metric  parameter  space  are  similar  to  those  for  the  situation 
where  the  parameter  space  is  a  general  complex  Banach  space.  When  deriving  the  sensitivity 
theory  for  the  convex  metric  parameter  space  case,  the  directional  derivative  is  used  instead 
of  the  Frechet  derivative  with  respect  to  the  measures. 

Given  any  two  arbitrary  points  fi\,u  G  (A4i,  djvq))  we  define  the  directional  derivative 
Sg(t,  x,  v  —  Hi)  of  g  at  pi  in  the  direction  v  —  y\  to  be  the  value  of  the  limit 


lim 

€■ — >0 
e>0 


g(t,  x,  //I  +  e(v  -  Ml))  -  g{t,  x,  Hi) 
e 


6g(t,x,ni,v  -  g i), 


provided  this  limit  exists  in  X.  A  sensitivity  theory  for  a  convex  metric  parameter  space 
as  developed  in  [13]  utilizes  the  directional  derivative.  We  make  the  assumptions  that  the 
function  g(t,  x,  pi)  of  (4.16)  has  a  bounded  smooth  Frechet  derivative  gx(t ,  x,  y\)  with  respect 
to  x  and  also  has  a  continuous  directional  derivative  6g(t,x,p, \,v  —  p\)  with  respect  to  p,\ 
in  the  direction  of  (i/  —  /xj).  Then  the  directional  derivative  y(t)  =  5x(t,  x ,  \  v  —  y,\)  exists, 

with  j/:R+xArxA4i— *  X,  and  y  satisfies  the  equation 


y(t)  =  gx(t,x(t,t0,x0,Hi),iJ.i)y(t)  +  6g(t,x(t,to,x0,fii),ni-,v  -  yx), 
y(to)  =  0.  (4.17) 

The  theories  of  [13,  27]  are  sufficient  to  a  develop  computationally  useful  sensitivity 
framework  and  also  provide  a  foundation  for  the  function  space  generalized  sensitivity  and 
asymptotic  error  analysis  proposed  in  the  next  two  sections. 


4.4  Generalized  Sensitivity 

Generalized  sensitivity  functions  (GSFs)  were  introduced  recently  by  Thomaseth  and  Cobelli 
[59]  as  a  new  tool  in  identification  studies  to  analyze  the  distribution  of  the  information 
content  (with  respect  to  the  model  parameters)  of  the  output  variables  of  a  system  for  a 
given  set  of  observations.  More  precisely,  the  generalized  sensitivity  functions  describe  the 
sensitivity  of  the  parameter  estimates  obtained  from  an  inverse  problem  with  respect  to  the 
data  measurements  used  in  the  inverse  problem. 

For  a  scalar  observation  model  with  discrete  time  measurements  (i.e.,  when  M  =  1  and 
C  is  a  1  x  Af  array  in  (4.10)),  the  generalized  sensitivity  functions  (GSFs)  are  defined  as 

*  1 

=  E  T5?Tt[f_1  x  V*/(t<,l9o)]  •  Vof(tu6 0),  (4.18) 

i=l  0  '  *' 

where  {ti},  l  =  1, . . . ,  n,  is  the  time  distribution  where  the  measurements  are  taken, 

F  =  E  ^7rrV,/(tj,0o)V,/(tj,0o)T  (4.19) 

;=1  °  Hi/ 
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is  the  corresponding  Fisher  information  matrix,  and  90  is  the  usual  “true”  parameter  value 
tacitly  assumed  to  exist  in  standard  statistical  theories.  The  symbol  represents  element- 
by-element  vector  multiplication.  For  the  motivation  and  details  which  led  to  the  definition 
above,  the  reader  should  consult  [59]  and  [34].  The  Fisher  information  matrix  measures  the 
information  content  of  the  data  corresponding  to  the  model  parameters.  In  (4.18)  we  see  that 
this  information  is  transferred  to  the  GSFs,  making  them  appropriate  tools  to  indicate  the 
relevance  of  the  measurements  to  the  estimates  obtained  in  parameter  estimation  problems. 

We  note  that  the  generalized  sensitivity  functions  (4.18)  are  vector- valued  functions  hav¬ 
ing  the  same  dimension  as  the  parameter  9  being  estimated.  The  fc-th  component  of  the 
vector  function  gs  represents  the  generalized  sensitivity  function  with  respect  to  9k-  The 
GSFs  are  defined  only  at  the  discrete  time  points  {tj,j  =  1, . . . ,  n}  and  they  are  cumulative 
functions,  at  each  time  point  f/  taking  into  account  only  the  contributions  to  estimates  of 
the  measurements  up  to  thus  indicating  the  influence  of  measurements  up  to  that  time 
on  the  parameter  estimates  obtained. 

From  (4.18)  it  is  readily  observed  that  all  the  components  of  gs  are  one  at  the  end  of 
the  total  time  interval,  i.e.,  gs (tM)  =  1.  If  one  defines  gs(t)  =  0  for  t  <  ti  (gs  is  zero 
when  no  measurement  is  collected),  then  each  component  of  gs  varies  from  0  to  1  during 
the  experiment.  As  developed  in  [59],  the  time  subinterval  during  which  this  transition 
has  the  sharpest  increase  corresponds  to  measurements  which  provide  the  most  information 
on  possible  variations  in  the  corresponding  estimated  model  parameters.  The  amount  of 
information  with  respect  to  a  parameter  9k  is  directly  related  to  the  rate  of  change  of  the 
corresponding  GSF;  thus  sharp  increases  in  gsfc  indicate  a  high  rate  in  additional  information 
about  9k  being  provided  by  the  new  measurements  in  that  time  period. 

For  finite  dimensional  systems  ( N  <  oo)  with  a  finite  number  ( p  <  oo)  of  parameters,  the 
numerical  implementation  of  the  generalized  sensitivity  functions  (4.18)  is  straightforward, 
since  the  gradient  of  /  with  respect  to  9  (or  xo)  is  simply  the  Jacobian  of  x  with  respect  to 
9  (or  Xo)  multiplied  by  C.  These  Jacobian  matrices  can  be  obtained  by  numerically  solving 
the  sensitivity  ODE  system  (4.13)  or  (4.14)  coupled  with  the  system  (4.9). 

We  have  effectively  used  GSFs  in  several  problems  [5,  12]  with  finite  dimensional  pa¬ 
rameters  in  ordinary  differential  equation  systems.  The  primary  contribution  of  generalized 
sensitivity  functions  is  that  they  can  indicate  regions  of  high  information  content  where, 
if  additional  data  points  are  taken,  one  can  generally  improve  the  existing  parameter  esti¬ 
mates.  This  has  rather  obvious  applications  to  experimental  design.  Moreover,  by  visually 
investigating  the  dynamics  of  GSF  curves,  one  can  potentially  identify  subsets  of  parameters 
which  are  highly  correlated.  While  there  are  still  some  heuristics  underlying  the  development 
of  GSFs  (the  presentation  in  the  appendix  of  [34]  offers  the  most  rigorous  presentation  to 
date),  it  is  clear  that  some  version  of  the  GSF  theory  should  become  a  valuable  tool  for  sci¬ 
entific  and  engineering  investigators  needing  to  estimate  parameters  in  dynamical  systems. 
Although  the  GSF  theory  provides  useful  information  in  parameter  estimation,  the  most 
insight  can  be  gained  when  the  GSFs  are  used  in  conjunction  with  traditional  sensitivity 
functions. 
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4.5  Asymptotic  Error  Estimates 

The  traditional  sensitivity  functions  that  we  have  introduced  above  also  play  a  fundamental 
role  in  quantifying  uncertainty  in  estimates  obtained  in  inverse  problems,  specifically  in  the 
computation  of  standard  errors  and  confidence  intervals.  We  consider  the  ordinary  least 
square  (OLS)  inverse  problem  by  seeking  9  that  minimize,  subject  to  (4.9),  the  cost  criterion 

m-'Eto-w)?,  (4.20) 

j= i 

where  {yj}  denotes  the  experimental  data. 

Once  the  optimal  6  are  found  using  an  optimization  algorithm,  standard  errors  and 
confidence  intervals  are  computed  by  using  the  asymptotic  theory  which  we  proceed  to 
outline.  Assume  n  scalar  measurements  (i.e.,  M  =  1)  are  represented  by  the  statistical 
model  (4.10)  or,  equivalently 

Yi  =  fj{e 0)+ej,  j  =  1, 2, ...  n,  (4.21) 

where  fj(9o)  =  f{tj,90)  is  the  model  for  the  observations  in  terms  of  the  state  variables  and 
0q  €  Mp  is  a  “set”  of  theoretical  “true”  parameter  values  (assumed  to  exist  in  a  standard 
statistical  approach).  Assume  further  that  the  errors  Cj,  j  =  1,2,  ...,n,  in  the  statistical 
model  of  the  observation  or  measurement  process  (4.21)  are  independent  identically  dis¬ 
tributed  ( i.i.d .)  random  variables  with  mean  E[tj]  —  0  and  constant  variance  var[tj ]  =  <7q 
where,  of  course  <7%,  is  unknown.  (The  constant  variance  assumption  can  be  validated  by 
use  of  standard  residual  plots  with  the  data  used  in  the  inverse  problems.)  It  follows  that 
the  observations  Yj  are  i.i.d.  with  mean  E\Yj\  —  fj(9o)  and  variance  var\Yj\  = 

Using  the  data  {yj}  for  the  observation  process  {V}}  with  the  model,  J(6)  is  optimized 
by  finding  the  OLS  estimator  6  in  (4.20).  Note  that  because  Yj  is  a  random  variable,  the 
estimator  <?ols  is  also  a  random  variable  with  a  distribution  called  the  sampling  distribution. 
Knowledge  of  this  sampling  distribution  provides  uncertainty  information  (e.g.,  standard 
errors)  for  the  numerical  values  of  6  obtained  using  a  specific  data  set  { y } }  (i.e.,  a  realization 
of  {V}})  when  minimizing  J(9). 

Under  reasonable  assumptions  on  smoothness  and  regularity  (the  smoothness  require¬ 
ments  for  model  solutions  are  readily  verified  using  continuous  dependence  results  for  or¬ 
dinary  differential  equations  in  many  examples;  the  regularity  requirements  involve,  among 
others,  conditions  on  how  the  observations  are  taken  as  sample  size  increases,  i.e.,  as  n  — ♦  oo), 
the  standard  nonlinear  regression  approximation  theory  ([39],  [45],  [55],  and  Chapter  12  of 
[58])  for  asymptotic  (as  n  — >  oo)  distributions  can  be  invoked.  This  theory  yields  that 
the  sampling  distribution  for  the  estimator  6(Y)  =  <?ols(}0>  where  Y  =  {Vj}"=1,  is  ap¬ 
proximately  a  p-multivariate  Gaussian  with  mean  E[9(Y)]  «  6q  and  covariance  matrix 
cov[9(Y)\  «  S0  =  cro[xT(0o)x(0o)]-1-  Here  x(0)  =  Fe(9)  is  the  n  x  p  sensitivity  matrix 
with  elements 

Xjk(0)  =  and  Fo(0)  =  . . . ,  fne(9))T. 

OVk 
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That  is,  for  n  large,  the  sampling  distribution  approximately  satisfies 


e0LS(Y)  ^  um(e0,a20[xT(e0)x(e0)}-1)  :=Mm(e0^0).  (4.22) 


The  elements  of  the  matrix  x  =  {Xjk)  can  be  estimated  using  the  forward  difference 


Xlk(e)  -~mk - N - ■ 

where  hk  is  a  p- vector  with  nonzero  entry  in  only  the  kth  component,  or  using  sensitivity 
equations  (see  [27]  and  the  references  therein).  Here  we  consider  the  sensitivity  equation 
approach  and  since  90,  a0  are  not  known,  we  approximate  them  in  S0  =  cro[xT(^o)x(^o)]_1- 
Following  standard  practice,  S0  is  approximated  by 

S0  »  S(0)  =  d2[xT(0)x(0)]_1 

where  9  is  the  parameter  estimate  obtained  from  minimizing  (4.20)  and  x($)  =  | £(0).  From 
the  outputs  defined  in  (4.9),  it  suffices  to  compute  the  sensitivities  ||  by  solving  the  N  x  p 
matrix  linear  variational  differential  equation 


d  / dx\  dgdx  dg 
di[d9j  =  frcd9  +  09' 


(4.23) 


which  is  the  same  as  the  sensitivity  equation  (4.13).  The  matrix  coefficient  and  the  forcing 
function  in  this  equation  are  evaluated  along  solutions  of  the  system  equation  (4.9).  The 
approximation  a2  to  a\  is  given  by 


p  j= 1 


Standard  errors  to  be  used  in  confidence  interval  calculations  are  thus  given  by  SEk{9 )  = 
y/Zkk(6),  k  =  l,2,...,p  (see  [36]). 

Finally,  in  order  to  compute  the  confidence  intervals  (at  the  100(1  —  c)%  level)  for  the 
estimated  parameters  <5C  and  ac,  we  define  the  confidence  level  parameters  associated  with 
the  estimated  parameters  so  that 

P(9k  -  tc/2SEk(9 )  <9k<9k  +  tc/2SEk{9 ))  =  1  -  c,  (4.24) 

where  c  €  [0, 1],  and  tc/2  €  M+.  For  a  given  c  value  (small,  say  c  =  .05  for  95%  confidence 
intervals),  the  critical  value  tc/2  is  computed  from  the  Student’s  t  distribution  tn~p  with  n  —  p 
degrees  of  freedom.  The  value  of  tc/2  is  determined  by  P(T  >  tc/2)  =  c/2  where  T  ~  tn~p. 

Again,  the  theory  outlined  above  is  for  finite  dimensional  systems  with  finite  dimensional 
parameters  9.  The  challenge  is  to  carry  out  an  asymptotic  analysis  in  infinite  dimensional 
function/measure  spaces  that  parallels  that  for  finite  dimensions  give  in  Chapter  12  of  [58]. 
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We  have  begun  using  the  approximation  ideas  of  [7,  13,  27,  30]  along  with  the  finite  dimen¬ 
sional  asymptotic  error  analysis  theory  outlined  above  as  applied  to  finite  nodal  values  to 
obtain  confidence  intervals  for  the  nodes,  leading  to  piecewise  smooth  “confidence  bands” 
or  confidence  functions  for  the  functions/distributions  being  estimated.  The  convergence 
theories  for  these  approximations  will  be  used  to  attempt  to  obtain  a  theory  for  limiting 
confidence  bands  for  the  infinite  dimensional  parameters  being  estimated.  We  have  some 
early  computational  results  providing  promise  that  this  research  could  lead  to  significant 
new  and  useful  methodologies.  In  [11]  we  investigated  stability  and  convergence  of  the  both 
delta  measure  and  spline  based  schemes  for  approximation  of  probability  distributions  (in  a 
Prohorov  framework  which  guarantees  convergence  of  the  distributions  but  not  the  densities 
in  all  cases-note  the  left  graphs  in  the  figures  below).  This  early  work  has  been  followed 
with  initial  efforts  at  obtaining  computationally  the  confidence  bands.  Sample  findings  are 
depicted  in  Figures  4  7  for  estimating  a  known  distribution  using  either  delta  measure  or 
spline  approximations  along  with  the  finite  dimensional  asymptotic  theory  to  compute  nodal 
confidence  intervals  for  the  densities.  These  nodal  confidence  intervals  are  then  interpolated 
to  produce  confidence  bands  for  both  the  densities  and  distributions.  These  initial  calcula¬ 
tions  suggest  existence  of  an  underlying  framework  for  asymptotic  confidence  bands  in  the 
case  of  infinite  dimensional  states  and  parameters. 

Known  and  Estimated  Probability  Densities  with  Truncated  Confidence  Intervals  (M«16)  Known  and  Estimated  Probability  Distributions  with  Confidence  Bands  (M*16) 


Figure  4:  Estimates  (red  lines)  of  probability  densities  and  probability  distributions  with 
confidence  intervals  and  bands  (dashed  lines)  given  a  true  gaussian  distribution  (blue  line) 
using  delta  measure  approximations  (DEL(16))  with  20%  absolute  error  in  data. 
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Known  and  Estimated  Probability  Densities  with  Truncated  Confidence  Intervals  (M*12)  Known  and  Estimated  Probability  Distnbutions  with  Confidence  Bands  (M«12) 


Figure  5:  Estimates  (red  lines)  of  probability  densities  and  probability  distributions  with 
confidence  intervals  and  bands  (dashed  lines)  given  a  true  gaussian  distribution  (blue  lines) 
using  spline  approximations  (SPL(12,128))  with  20%  absolute  error  in  data. 


Known  and  Estimated  Probability  Densities  with  Truncated  Confidence  Intervals  (M=28) 


Figure  6:  Estimates  (red  lines)  of  probability  densities  and  probability  distributions  with 
confidence  intervals  and  bands  (dashed  lines)  given  a  true  bi-gaussian  distribution  (blue 
lines)  using  delta  measure  approximations  (DEL(28))  with  20%  absolute  error  in  data. 
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Known  and  Estimated  Probability  Densities  with  Truncated  Confidence  Intervals  (M*=12)  Known  and  Estimated  Probability  Distributions  with  Confidence  Bands  (M-12) 


Figure  7:  Estimates  (red  lines)  of  probability  densities  and  probability  distributions  with 
confidence  intervals  and  bands  (dashed  lines)  given  a  true  bi-gaussian  distribution  (blue 
lines)  using  spline  approximations  (SPL(12,128))  with  20%  absolute  error  in  data. 
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4.6  Nodal  Networks  with  Uncertainty 

Networks  with  uncertainty  are  ubiquitous  in  science  and  engineering,  but  especially  in  com¬ 
plex  nonlinear  systems  of  interest  to  DOD.  In  particular,  nodal  networks  with  uncertainty  are 
fundamental  to  the  understanding  of  problems  involving  dynamic  resource  management  and 
allocation.  Examples  range  from  service  scheduling  and  material  supply  chains  to  produc¬ 
tion  systems  and  modern  high  speed  communication  networks  where  information  is  routed 
between  nodes  in  a  distributed  information  system  with  a  global  grid. 

We  summarize  one  approach  in  an  interdisciplinary  group  effort.  In  [5],  an  approach  to 
modeling  the  impact  of  disturbances  in  an  agricultural  production  network  is  presented.  A 
stochastic  model  and  its  approximate  deterministic  model  for  averages  over  sample  paths 
of  the  stochastic  system  are  developed.  Simulations,  sensitivity  and  generalized  sensitivity 
analyses  were  given.  Finally,  it  was  shown  how  diseases  (our  example  is  modeled  on  the 
devastating  Foot  and  Mouth  Disease  (FMD)  that  hit  the  United  Kingdom  (UK)  in  2001) 
may  be  introduced  into  the  network  and  corresponding  simulations  are  discussed. 

The  current  production  methods  for  livestock  follow  the  just-in-time  philosophy  of  many 
manufacturing  industries  as  well  as  supply  chains.  Feedstock  and  animals  are  grown  in 
different  areas.  Animals  are  moved  from  one  farm  to  another,  depending  on  their  age. 
Unfortunately,  shocks  propagate  rapidly  through  such  systems;  an  interruption  to  the  feed 
supply  has  a  much  larger  impact  when  farms  have  minimal  surplus  supplies  in-stock  than 
when  they  have  large  reserves.  The  just-in-time  movement  of  animals  between  farms  serves 
as  another  vulnerability.  Stopping  movement  of  animals  to  and  from  a  farm  with  animals 
infected  by  a  disease  will  have  effects  that  quickly  spread  through  the  system.  Nurseries 
supplying  the  farm  will  have  nowhere  to  send  their  animals  as  they  grow  up.  Finishers  and 
slaughterhouses  supplied  by  the  farm  will  have  their  supply  interrupted. 

In  [5]  we  combined  initial  statistical  and  mathematical  modeling  ideas  to  address  the 
above  issues,  using  the  North  Carolina  swine  industry’s  potential  response  to  FMD  as  an 
example.  We  focused  our  attention  on  the  North  Carolina  swine  industry  because  it  is  both 
the  second  largest  swine  industry  in  the  United  States,  and  is  local  to  us.  Our  goal  was  to 
develop  models  that  could  be  used  to  investigate  how  small  perturbations  to  the  agricultural 
supply  system  would  affect  its  overall  performance.  A  hurricane  that  throttles  inter- farm 
transportation  for  a  short  duration,  or  a  disease  outbreak  that  spreads  through  distribution 
channels  are  example  causes  of  the  perturbations  of  interest.  In  the  former  case,  the  just-in- 
time  delivery  systems  may  not  provide  enough  slack  to  absorb  the  shock.  In  the  latter  case, 
strategies  that  involve  destruction  of  all  livestock  in  an  infected  branch  of  the  system  may 
be  overly  harsh;  a  more  moderate  response  may  be  as  effective  without  the  high  toll  on  the 
infrastructure. 

We  modeled  a  simplified  swine  production  network  in  North  Carolina  containing  four 
levels  of  production  nodes:  growers/sows  ( Node  1),  nurseries  ( Node  2),  finishers  ( Node  3), 
and  processing  plants/slaughter  houses  ( Node  4).  At  grower  or  sow  farms  ( Node  1),  the  new 
piglets  are  born  and  typically  weaned  three  weeks  after  birth.  The  three- week  old  piglets  are 
then  moved  to  the  nursery  farms  ( Node  2)  to  mature  for  another  seven  weeks.  They  are  then 
transferred  to  the  finisher  farms  ( Node  3),  where  they  grow  to  full  market  size,  which  takes 
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about  twenty  weeks.  Once  they  reach  market  weight,  the  matured  pigs  are  moved  to  the 
processors  (slaughter  houses)  ( Node  4).  Pork  products  then  continue  through  wholesalers  to 
consumers.  There  are  also  several  inputs  to  the  system  which  we  will  not  consider,  such  as 
food,  typically  corn  grown  in  the  midwest.  There  are  several  types  of  breeding  farms  where 
pure-bred  stock  are  raised;  these  are  typically  crossed  to  produce  hybrid  strains  for  pork 
production. 

In  our  efforts  in  [5],  we  formulated  a  nonlinear  stochastic  model  for  our  agricultural 
network  and  showed  how  it  can  be  converted  to  an  equivalent  (in  a  sense  made  precise 
below)  deterministic  differential  equation  model.  This  deterministic  model  readily  lends 
itself  to  simulations  and  as  well  as  to  sensitivity  analysis  techniques  such  as  those  outlined 
above.  In  particular,  we  used  it  in  simulations  to  mimic  the  presence  of  an  infectious  disease 
in  the  network. 

In  addition  to  the  development  of  models  for  a  typical  production  network  permitting 
perturbations,  a  significant  contribution  of  the  effort  in  [5]  is  the  demonstration  of  stochastic, 
mathematical  and  computational  methodology  that  is  available  to  domain  scientists,  statisti¬ 
cians  and  applied  mathematicians  working  in  a  concerted  team  effort  on  complex  problems 
of  the  type  exemplified  here. 

4.7  Modeling 

We  considered  in  [5]  a  simplified  swine  production  network  with  four  aggregated  nodes:  sows 
(Node  1),  nurseries  (Node  2),  finishers  (Node  3),  and  slaughter  houses  (Node  4).  Our  goal 
was  to  study  the  effects  of  perturbations  within  the  network.  This  can  be  done  either  by 
affecting  the  nodes  or  the  transitions  between  nodes  directly  or  indirectly.  For  instance,  a 
problem  with  the  breeding  farms  would  result  in  a  reduction  of  sows  available  for  producing 
new  piglets.  This  would  result  in  a  reduced  rate  of  transition  from  Node  1  to  Node  2, 
since  we  could  not  grow  as  many  piglets.  We  could  then  track  the  effect  of  this  through  our 
network. 

Although  unavoidable  in  actual  production  processes,  we  assumed  in  our  example  that 
there  are  no  net  losses  in  the  network  (i.e.,  the  total  number  of  pigs  in  the  network  remains 
constant)  and  that  the  only  deaths  occur  at  the  slaughter  houses.  Thus  we  assumed  that  the 
number  of  processed  pigs  per  day  at  the  slaughter  houses  is  equal  to  the  number  of  newborn 
piglets  per  day  at  the  growers.  We  can  model  reduced  birth-rates  by  reducing  the  rate  at 
which  piglets  move  to  the  nurseries.  This  led  us  to  deal  with  a  closed  network.  We  note 
that  this  approximation  is  realistic  when  the  network  is  efficient  and  operates  at  or  near  full 
capacity  (i.e,  when  the  number  of  animals  removed  from  the  chain  are  immediately  replaced 
by  new  production/growth,  avoiding  significant  idle  times).  Our  closed  network  model  for 
the  swine  production  is  summarized  schematically  in  Figure  8. 

Each  node  with  corresponding  population  number  Ni,  i  =  1, . . . ,  4,  in  Figure  8  represents 
an  aggregation  of  all  the  production  units  corresponding  to  that  level  in  the  production 
network.  Given  a  specific  production  network,  any  of  the  four  levels  of  the  chain  may  be 
broken  into  its  constituent  units  (e.g.,  farms),  and  analyzed  in  detail  as  a  separate  sub- 
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Figure  8:  Aggregated  agricultural  network  model. 


network.  The  directed  edges  between  the  nodes  represent  the  movement  of  the  pigs  through 
the  network.  The  rate  is  determined  by  the  pigs’  residence  time,  the  number  of  pigs  at  each 
node,  and  the  capacity  constraints  at  the  corresponding  nodes.  Let  L,  denote  the  capacity 
constraint  at  node  z,  i  =  1, ...  ,4.  Since  we  have  a  closed  network,  it  is  assumed  that  there 
is  no  capacity  constraint  at  Node  1  and  therefore  we  take  L\  =  oo.  We  also  define  Sm  to 
be  the  maximum  exit  rate  at  Node  4,  i.e. ,  the  maximum  killing  capacity  at  the  slaughter 
house.  The  residence  times  at  each  node,  together  with  the  capacity  constraints  and  the 
slaughterhouse  killing  capacity,  are  based  on  very  rough  estimates  of  swine  production  in 
North  Carolina  and  these  were  used  to  obtain  parameters  for  the  simulations  reported  on 
here  and  in  [5]. 


4.8  Stochastic  and  Deterministic  Models 


We  model  the  evolution  of  the  food  production  network  shown  in  Figure  8  as  a  continuous 
time  discrete  state  density  dependent  jump  Markov  Chain  ( MC j  [3,  42]  with  a  discrete  state 
space  embedded  in  an  K4  non-negative  integer  lattice  C.  The  state  of  this  MC  at  time  t  is 
denoted  by  X(f)  =  (Xi(f), . . . ,  X4(t)),  where  Xi(t)  is  the  number  of  pigs  at  node  i  at  time 
t,i  =  1,...,4. 

The  rates  of  transition  of  X(£)  are  nonlinear  functions  A  i  :  C  —*  [0,  oo)  for  i  =  1, . . . ,  4, 
and  for  x  €  C  are  given  by: 


Ai(x)  :=  qi(xx  -  l,x2  +  1, ar3, a:4) 
A2(x)  :=  q2(x i,x2  -  l,x3  + l,x4) 
A3(x)  :=  q3(xi,x2,x3  -  l,x4  +  1) 
A4(x)  :=  q4{x  1  +  1,  x2,  x3,  x4  -  1) 


fciXj  (L2  —  x2)+ 
k2  x2  (L3  -  x3)+ 
k3  x3  {La  -  x4)+ 

k4  min(x4,  Sm )  (4-25) 


where  fcj,i  =  1,...,4,  is  proportional  to  the  service  rate  at  node  i;  Lj,i  =  2,3,4,  is  the 
buffer  size  (capacity  constraint)  at  node  i  and  Sm  is  the  slaughter  capacity  at  node  4  as 
discussed  above.  For  any  real  z,  the  symbol  (z)+  is  defined  as  the  non-negative  part  of  z, 
i.e.,  (z)+  =  max(z, 0).  Then  q\{xi  —  l,x2  +  l,x3,x4)  is  given  by 


<7i(xi  -  l,x2  +  l,x3,x4)  = 

,.  Pr[X(f  +  h)  =  (xi  -  l,x2  +  1, x3, x4)|X(f)  =  (xi, 
lim  — 1 - - - 

h — ►0-j-  h 


X2,X3,XA)) 
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The  other  g*  are  given  similarly. 

The  simple  model  for  transition  rates  given  above  is  formulated  under  the  following 
assumptions  and  hypotheses.  First  it  is  assumed  that  the  transportation  rates  qt,i  —  1,2,3, 
are  proportional  to  x,  (L,+1  —  Xj+i)+,  the  product  of  the  number  of  animals  available  and  the 
available  capacity  at  the  next  node.  If  no  capacity  is  available,  the  rate  is  taken  as  zero.  The 
rate  at  the  slaughter  house  (Node  4)  is  the  maximum  Sm  if  a  sufficient  number  of  animals  is 
available;  otherwise  all  animals  present  at  the  node  are  slaughtered  on  that  day.  Finally,  it 
is  assumed  that  the  network  is  at  or  near  steady-state  and  maximum  efficiency  in  that  the 
slaughter  rate  at  Node  4  is  the  same  as  the  input  at  Node  1  (this  is  represented  schematically 
in  Figure  8  by  the  arrow  from  Node  4  to  Node  1).  This  results  in  the  rate  dynamics  (4.32) 
below  with  the  output  rate  at  Node  4  the  same  as  the  input  rate  at  Node  1. 

We  remark  that  the  product  nonlinearities  x,  (Li+ 1  — xi+i)+  of  (4.25)  where  transportation 
occurs  more  rapidly  the  further  the  node  level  is  from  capacity  (i.e.,  the  system  reacts  more 
rapidly  to  larger  perturbations  from  capacity)  are  only  one  possible  form  for  these  terms. 
One  could  also  reasonably  argue  for  alternative  terms  of  the  form  x*  Xi+i  where  xl+i  is  the 
characteristic  function  for  the  set  {(L,+1  —  xi+1)  >  0}  so  that  the  transportation  rate  from 
a  node  depends  only  on  the  number  available  at  that  node  so  long  as  capacity  at  the  next 
node  has  not  been  reached.  We  note  that  in  this  case  the  sensitivity  analyses  using  the  ideas 
discussed  above  are  more  difficult  due  to  a  lack  of  continuity  of  the  dynamics  in  the  system 
equations. 

Let  Ri(t)  i  =  1, . . . ,  4,  denote  the  number  of  times  that  the  ith  transition  occurs  by  time 
t.  Then  Ri  is  a  counting  process  with  intensity  A j(X(t)),  and  the  corresponding  stochastic 
process  can  be  defined  by 

Ri(t)  =  Yi{[  MX(s))ds),  i  =  1 . . . ,  4,  (4.26) 

where  the  Yt  are  independent  unit  Poisson  processes.  That  is,  sample  paths  r,(f)  of  Rl(t) 
are  given  in  terms  of  sample  paths  x(i)  of  X(t)  by 

n(t)  =  Yi(J  \i(x(s))ds),  *  =  1  —  ,4.  (4.27) 

We  write  Ri  in  this  form  to  illustrate  that  A*  is  a  rate  of  the  corresponding  counting  process. 

Our  stochastic  model  then  can  be  written  as 

*i(0  =  X^O)  -  Ri(t)  +  R<(t) 

X2(t)  =  *2(0)  +  Ri(t)-R2(t) 

Xz(t)  =  X3(0)  -I-  R2(t)  —  Rs(t ) 

X4(t )  =  X4(0)  +  R3(t)  -  R4(t).  (4.28) 

The  above  system  typically  cannot  be  solved  for  a  stationary  distribution  and  an  empirical 
approach  based  on  the  so-called  Gillespie  algorithm  [47]  can  be  used  to  investigate  the 
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long  term  behavior  of  the  system  (see  [5]).  The  approximate  large  population  behavior  of 
an  appropriately  scaled  system  may  be  also  analyzed  via  macroscopic  deterministic  rate 
equations  as  we  shall  explain  next  (the  original  theory  is  due  to  Kurtz  and  is  discussed  in 
[42]  and  the  references  therein). 

Let  N  be  the  total  network  or  population  size.  If  N  is  known  we  may  consider  the 
animal  units  per  system  size  or  the  units  concentration  in  the  stochastic  process  CN(t)  = 
X.(t)/N  with  sample  paths  cN(t).  For  large  systems  this  approach  leads  to  a  deterministic 
approximation  (obtained  as  solutions  to  the  system  rate  equation  defined  below)  to  the 
stochastic  equation  (4.28),  in  terms  of  c(£),  the  large  sample  size  average  over  sample  paths 
or  trajectories  c N(t)  of  CN(t). 

We  rescale  the  rate  constants  k,  Li  and  Sm  as  follows: 

k4  =  k i,  Ki  =  Nki,  i  =  1, 2  or  3, 

=  Sm/N,  h  =  Lt/N.  (4.29) 

According  to  Equation  (4.25),  this  rescaling  implies  that 


A i(x)  =  KiXi(Li+ 1  -  Xi)+/N  =  NKiC?(li+1  -  cf)+  i  =  1,2,3, 

and 

A4(x)  =  k4  min(x4,S'm)  =  Nk4  min(c^,sm). 

Recall  that  for  large  N  the  Strong  Law  of  Large  Numbers  (SLLN)  for  the  Poisson  Process 
Y  implies  Y ( Nu)/N  ~  u  [48].  One  can  use  this  fact,  along  with  the  rescaling  of  the  constants 
as  given  above,  to  argue  that  sample  paths  rf  t)  for  the  counting  process  (4.26)  defined  in 
terms  of  the  sample  paths  x(£)  or  cN(t)  =  x.(t)/N  may  be  approximated  for  large  N  in  terms 
of  the  deterministic  variables  c (t),  the  averages  over  sample  paths  or  trajectories  c N(t)  of 
C N{t),  by 


r\N\t)  =  =  JjYi(f0  Xi(x(s))ds) 

=  ±Yi(Nj\ic»(s)(li+1-cUs))+  ds) 

KiCi(s)(li+1  -  d+i(s))+ds  fori  =  1,2, 3,  (4.30) 

and  similarly 

1  f* 

r{N) (t)  = —r4(t)  &  J  k4  min(c4(s),  sm)  ds. 

For  a  full  and  rigorous  discussion  of  this  “approximation  in  mean”,  see  Chapters  6.4  and  11  of 
[42]  and  Chapter  5  of  [3].  (We  note  that  these  modeling  ideas  are  also  useful  in  examples  for 
logistic  growth,  epidemics  (SIR),  chemical  reaction  networks  and  general  “density  dependent 
processes”,  among  others.)  The  averages  c(t)  satisfy  a  system  of  deterministic  ordinary 
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differential  equations  which  can  be  heuristically  derived  by  beginning  with  the  system  (4.28). 
Upon  dividing  both  sides  of  each  equation  by  N  and  applying  the  above,  we  obtain  the  rate 
equations ,  i.e.,  the  system  of  integral  equations  approximating  via  the  SLLN  the  original 
stochastic  system,  as  follows: 


cT(t) 


cf(t) 


ci(0) 

r[N\t )  +r[N) 

w 

r  t 

cj(0) 

— 

rt 

/  KiCi(s){12 
Jo 

-c2(s))+ 

ds 

+ 

/  /t4  min(c4(s), 
)o 

•^m)  ds 

c2(0) 

+  ' 

-TW  -  r<"» 

w 

c2(  0) 

— 

rt 

/  K2C2(s)(/3 
Jo 

-  c3(s))  + 

ds 

+ 

rt 

/  «lCi(s)(/2- 
Jo 

c2{s))+ds 

ca(0) 

+ 

r no  -  -r 

w 

c3(0) 

— 

j  k3c3(s)(14 

-  c4(s))+ 

ds 

+ 

rt 

/  K2C2(s)(/3  — 
Jo 

c3{s))+ds 

c4(0) 

+ 

rfrm-rr 

(t) 

rt 

c4(0) 

+ 

rt 

/  K3C3(s)(/4 
Jo 

-  C4(s))  + 

ds 

/  /t4min(c4(s), 
Jo 

i  ®m)  ds. 

Upon  approximating  the  c?(t)  on  the  left  above  by  the  c,(£)  and  differentiating  the 
resulting  equations,  we  find  that  the  integral  equation  system  is  equivalent  to  a  system  of 
ordinary  differential  equations  for  c(t)  €  R4  given  by 


/CiCi(i)(Z2  -  c2(t))+  +  K4min(c4(£),  sm) 

ft2C2(£)(^3  —  c3(<))+  +  «lCi(f)(/2  —  C2(f))  + 

«3C3(£)(/4  -  C4(£))+  +  «2c2(£)(/3  -  c3(t))+ 

«4min (c4(£),  sm)  +  k3c3(£)(/4  -  c4(£))+  (4.32) 

with  the  initial  conditions  c(0)  =  Co-  As  we  demonstrate  next,  solutions  of  these  equations 
yield  quite  good  approximations  to  the  sample  paths  of  the  stochastic  system. 

A  sample  of  five  realizations  with  N  =  3, 165,000  is  plotted  in  Figure  9(top).  Note  that 
the  realizations  exhibit  very  little  visible  differences.  However,  when  one  carries  out  the 
simulations  for  a  smaller  system  (N  =  3, 165  pigs  with  the  parameters  scaled  accordingly), 
the  variations  are  readily  visible  as  can  be  seen  in  Figure  9  (bottom).  We  also  remark  that 
these  simulations  offer  graphic  depictions  of  the  approximation  theory  discussed  above  where 
in  the  case  of  very  large  N  one  can  cannot  distinguish  between  the  stochastic  simulations  and 
the  corresponding  deterministic  simulations  for  the  sample  path  averages.  This  comparison 
is  given  in  Figure  10  where  we  have  depicted  (bottom)  the  solution  of  the  concentrations 
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Figure  9:  Simulation  of  the  stochastic  model  (4.28)  of  the  agriculture  production  chain  via 
the  Gillespie  algorithm  for  N=3, 165,000  (top)  vs.  N=3,165  (bottom) 
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in  system  (4.32)  as  the  rescaled  quantities  Ni(t)  =  Nci(t)  vs.  the  corresponding  stochastic 
simulation  (top).  As  expected,  we  find  that  the  stochastic  and  deterministic  computations 
provide  similar  numerical  results  with  the  realizations  fluctuating  about  the  solution  of  the 
deterministic  system  for  the  averages  as  predicted  by  the  theory. 
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Figure  10:  Simulation  with  the  stochastic  model  (top)  vs.  solution  of  the  deterministic  (rate 
equation)  system  (bottom)  for  t  €  [0, 100]  with  N=3, 165,000. 
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4.9  Summary  Comments  on  Nodal  Networks 

In  [5]  we  demonstrated  a  methodological  approach  to  the  investigation  of  production  net¬ 
works  and  their  vulnerability  to  disturbances  such  as  diseases.  The  stochastic  model  and 
the  resulting  approximate  deterministic  system  we  employed  were  shown  to  agree  well.  We 
also  carried  out  simulations  and  sensitivity  analyses.  We  used  the  deterministic  model  to 
show  how  to  determine  the  parameters  to  which  the  model  exhibits  the  most  sensitivity. 
This  is  more  readily  done  with  the  deterministic  approximation  of  the  stochastic  model,  but 
allows  one  to  draw  conclusions  about  sensitivity  with  respect  to  parameters  in  the  Markov 
process.  We  also  demonstrated  in  [5]  how  disease  can  be  introduced  and  the  resulting  net¬ 
work  vulnerabilities  analyzed.  The  deterministic  version  of  the  model  is  also  amenable  to 
the  use  of  inverse  problem  algorithms  with  the  data  to  obtain  parameter  estimates  along 
with  measures  of  associated  uncertainty  (e.g.,  standard  errors  [36,  58])  for  the  underlying 
transition  rates  A;*.  The  transmission  dynamics  of  infection  used  in  [5]  were  deterministic  in 
nature.  However,  it  is  well  known,  for  example,  that  random  effects  can  have  a  major  impact 
on  the  invasion  of  an  infection  into  a  population.  Finally,  the  model  as  developed  assumes 
instantaneous  transport  between  nodes.  If  infection  during  transport  is  an  important  factor 
(and  depending  on  the  disease  it  may  well  be)  then  the  structure  of  the  model  should  be 
modified  to  incorporate  positive  transport  times.  This  could  lead  to  more  interesting  (math¬ 
ematically)  and  more  difficult  dynamical  systems  with  time  delays  in  place  of  (4.32)  and 
the  corresponding  systems  in  [5].  All  of  these  issues  should  be  addressed  in  a  general  nodal 
network  model. 

The  randomness  seen  in  the  stochastic  network  model  originates  from  the  random  move¬ 
ment  of  discrete  individuals  from  node  to  node.  The  analysis  of  [5]  and  the  figures  above 
show  that,  due  to  an  averaging  effect,  these  random  effects  become  less  important  as  the 
system  size  N  increases.  Application  of  the  stochastic  transportation  model  to  describe  real- 
world  situations  should,  therefore,  account  for  the  size  of  the  groups  in  which  units  (pigs  in 
the  case  of  the  agriculture  production  model  of  [5])  are  transported  between  nodes.  If,  for 
example,  one  thousand  units  were  moved  at  a  time,  the  appropriate  notion  of  an  “individual” 
within  the  model  would  be  a  thousand  units.  Treating  each  group  of  a  thousand  individuals 
as  a  unit  would  lead  to  a  marked  increase  in  the  magnitude  of  stochastic  fluctuations  seen  at 
the  population  level.  Consequently,  even  though  the  system  size  in  simulations  might  be  on 
the  order  of  millions  of  units,  it  might  be  that  the  resulting  stochastic  fluctuations  in  a  more 
realistic  model  of  such  production  systems  are  closer  to  those  shown  in  Figure  9  (bottom) 
than  to  those  of  Figure  9  (top). 

The  approach  developed  in  [5]  and  outlined  in  this  paper  has  rather  obvious  potential 
for  application  to  a  wide  range  of  problems.  These  include  the  investigation  of  the  spread 
of  diseases  through  spatially  or  structurally  distributed  dynamic  populations  (e.g.,  avian 
flu  through  migrating  bird  populations,  contagious  infections  through  highly  mobile  and/or 
age-structured  human  and  animal  populations).  In  some  of  these  cases  the  natural  nodal 
structure  would  be  a  continuum,  requiring  stochastic  and  deterministic  models  with  a  con¬ 
tinuum  of  spatial/structural  heterogeneities,  leading  to  partial  differential  equation  systems. 
Such  applications  would  undoubtedly  motivate  the  development  of  interesting  new  stochastic 
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and  deterministic  mathematical  and  computational  methodologies. 

More  importantly,  the  approach  and  methodology  outlined  here  are  useful  in  general 
supply  networks  for  investigation  of  a  wide  range  of  perturbations  (either  accidental  or 
deliberate)-e.g.,  loss  of  capacity  at  a  given  node  such  as  a  factory  being  shut  down,  a  node 
being  rendered  inoperable  via  weather,  technical  (electrical/mechanical)  difficulties  in  com¬ 
munication/transport  systems,  etc. 
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